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INTRODUCTION 
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A Fortran  program  designated  as  USL  Program  No.  0293  has  been 
written  by  the  Digital  Computing  Branch  for  the  purpose  of  carrying  out 
a generalized  harmonic  analysis  of  a digitized  time  function.  This 
analysis  procedure,  as  well  as  the  program  itself,  will  be  described  in 
this  memorandum.  , , . . , , 

1 If  -■  ti-  ' -• 

ANALYSIS 


(r<7 
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A Fburier  transform  pair  relationship  exists  between  the  autocovariance 
function  and  the  power  spectral  density  of  a stationary  random  process 
x(t).  Using  this  association  a method  for  obtaining  an  estimate  of  the 
power  spectrum  of  the  process  by  transforming  the  autocovariance  function 
has  been  developed  by  Blackman  and  Tukey.1  Basically  this  analysis  consists 
of  first  calculating  the  normalized  autocovariance  function 


| f - 1 ^ ' r(T)  = x’(t)  x'(t+T)  (1) 

of  the  sample  function  x'(t)  and  modifying  this  function  using  the  Hamming 
lag  window2.  x'(t)  is  only  one  member  of  the  ensemble  of  functions  that 
make  up  the  random  process  x(t).  Consequently  it  is  only  possible  to 
obtain  estimates  of  the  power  spectral  density  of  the  random  process 
by  performing  a Fourier  transformation  of  this  modified  correlation  function. 


1 R.  B.  Blackman  and  J.  W.  Tukey,  "The  Measurement  of  Power  Spectra 
from  the  Point  of  View  of  Communications  Engineering",  1958 


2 Other  such  windows  are  discussed  but  the  Hamming  window  is  used  in 
the  program  described  by  this  memorandum. 
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Since  x'(t)  is  only  available  for  a finite  time  T,  r(T)  can  only 
be  calculated  for  /T/  =0,  Tm,  where  Tm  is  the  maximum  lag 

time,  and  ia  not  defined  outside  this  region.  The  Fourier  transform 
of  r(T)  is  given  by 


S(f) 


r ( T)  e-J“T  dT 


(2) 


In  order  to  carry  out  this  integration,  r(T)  must  be  defined  for  all 
T in  (-*0,  By  considering  r*  (T)  as  a product  of  r(T)  with  a 

weighting  function  or  lag  window  w(t)  where 


w(T) 


1 I T I < Tm 

0 otherwise 


(3) 


r'(T)  is  defined  for  all  T and  its  Fourier  transform  may  be  calculated. 
Since  multiplication  in  the  time  domain  corresponds  to  convolution  in 
the  frequency  domain,  the  resulting  estimated  power  spectrum  is  given  by 


P(f)  = S(f)  * W(f) 


OO 


where  S(f)  is  defined  by  eq.  (2)  and 


W(f) 


W(T)e-J“'TdT 


2 T„ 


(5) 


Vf(f)  has  very  undesirable  minor  lobes  which  contribute  errors  to  the 
power  spectral  density  estimates.  For  this  reason  other  weighting 
functions  have  been  examined  and  are  described  in  footnote  (1).  The 
particular  lag  window  being  used  in  Program  No.  0293  is  the  Hamming 
window  given  by 


h(T) 


.54  + .46 


0 


|t|  5.  T 


I T I > T 


(6) 
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A comparison  of  h(T)  and  w(T)  is  shown  in  Figure  1.  The  Hamming 
window  was  chosen  because  it  has  the  property  that  its  Fourier 
transform 

H(f)  = .54  W(f)  + .23W(f  + j.._)  +.23W(f-  1 1 (7) 


has  minor  lobes  that  are  no  more  than  2 per  cent  of  the  height  of 
the  major  lobe,  a vast  improvement  over  VJ(f ) which  has  side  lobes 
as  high  as  20  per  cent  of  the  major  lobe.  These  differences  are 
illustrated  in  Figure  2.  In  addition  to  this  important  property 
the  Hamming  window  is  easy  to  work  with  computationally. 

Once  having  selected  a weighting  Function,  power  spectral  density 
estimates  are  obtained  by  performing  a Fourier  transformation  on  the 
product  r(T)  * h(T).  Since  both  functions  are  even,  a cosine 
transformation  is  all  that  is  required. 

In  practice,  and  in  this  program,  continuous  data  is  not  available 
but  instead  discrete  values  of  the  sample  function  x'(t)  at  equal 
increments  At  are  available  for  analysis.  In  this  situation,  calculations 
are  simplified  by  the  following  approach.  First  form  mean  lagged 
products  r(Tp)  p = 0,  1,  ...  m,  (the  discrete  analog  of  r(T))  where 
Tm  = mAt . Next  perform  a finite  cosine  series  transform  using 

these  ra  + 1 correlation  coefficients  giving  raw  estimates  of  the 
power  spectral  densityEL(p).  These  estimates  are  taken  equally 
spaced  at  p = l/2Tm  = l/2mAt  . These  raw  estimates  correspond  to 
using  the  rectangular  weighting  function  w(t)  described  in  the  continuous 
case.  Just  as  a modification  was  required  in  the  continuous  case 
to  reduce  the  errors  imposed  by  the  large  minor  lobes  of  W(f),  so  too 
smoothed  estimates  of  the  power  spectral  density  must  be  calculated 
in  the  discrete  case.  The  only  major  difference  being  that  in  the 
previous  case  the  modification  was  carried  out  before  transforming 
in  order  to  avoid  convolution  in  the  frequency  domain.  With  discrete 
samples,  however,  computations  are  easily  carried  out  in  the  frequency 
domain.  Since  the  two  functions  being  convolved  are  just  impulses  at 
equal  spacings  convolution  results  in  another  set  of  impulses.  More 
specifically  for  each  raw  estimate  EL(p)  convolution  results  in 
smoothed  estimates  Sp  where 


S(p)  = .23  EL(p-l)  +.54  EL(p)  + .23  EL(p+l) 


(8) 
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RESOLUTION  AND  STABILITY 


Using  the  method  outlined,  estimates  of  the  power  spectral 
density  are  obtained  every  1 cps.  These  estimates  are  considered 

2Tm 

to  be  an  average  over  frequency  of  the  power  in  the  band  l/2Tm  cps. 

As  seen  in  equation  8,  and  in  Figure  2,  adjacent  estimates  are  not 
independent  so  the  resolution  is  taken  as  / 1 \ = 1 cps. 

^ 2Tm  ) Tm 


The  stability  of  these  estimates  is  related  to  the  number  of 
degrees  of  freedom  associated  with  each  estimate.  If  the  true  spectrum 
is  relatively  smooth,  degrees  of  freedom  may  be  approximated  by 

K 2T'  (9) 

Tm 


where  T'  = T - l/3  Tm. 

2 2 

Call  the  estimates  of  the  power  spectral  density  S and  let  o 
be  the  true  power  spectral  density.  The  random  variable  KS^/ <7  has 

a^  X2  distribution  with  K degrees  of  freedom.  Using  the  table  of 
x values,  it  is  possible  to  calculate  confidence  limits  for  <r  the 
true  power  spectral  density.  The  90$  confidence  limits  are  given  by 


< c 


where 


*=  2 2 2 2 2 
X-^  and  X2  are  such  that  P(X  > Xj_  ) = .95  and  P(X  > X2  ) • .05 

For  example,  if  K = 25  degrees  of  freedom,  then  90$  of  the  time  the 
true  estimate  will  fall  in  the  interval  given  by  ,66a  < <r  2 < 1.71s 

Since  the  length  of  the  sample  function  x'(t)  is  limited  to  T 
seconds  and  stability  is  dependent  upon  keeping  the  ratio  T/Tm  large, 
a compromise  is  always  necessary  between  stability  and  resolution. 

As  Tm  increases,  so  does  the  resolution,  but  the  stability  or 
confidence  in  the  estimates  of  the  power  spectrum  decreases.  For  a 
fixed  T,  & value  of  Tm  will  be  determined  for  each  analysis  to  give 
the  required  resolution  with  the  accuracy  needed,  however,  a rule  of 
thumb  is  that  Tm  should  never  be  greater  them  5 or  10  per  cent  of  T 
thus  always  allowing  at  least  20  degrees  of  freedom  for  each  estimate. 
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Suppose  it  is  desired  to  analyze  frequencies  from  0 xo  some 

f cps.  if  the  true  spectrum  is  zero  outside  f then  a sampling 
max  max 

rate  of  At  = 1 is  sufficient.  However,  if  the  data  is 

2 f 
max 

filtered  and  the  frequency  response  of  the  filter  drops  off  gradually 

outside  f then  perhaps  a sampling  rate  of  At  ■ 1 should  be 

max  « a 1 

j xmax 

used  to  guard  against  the  effect  of  aliasing  of  frequencies  above  fmax. 

Suppose  a resolution  of  1 cps  is  desired  and  additionally  the  accuracy 
of  50  degrees  of  freedom  is  sufficient,  then 


K = 2T1  = 50  and  1 = 1 cps 

Tm  Tm 


(11) 


Therefore  Tm  = 1 second  and  T*  = 25  seconds. and  T =* 25  seconds. 

This  means  that  if,  for  example,  f = 100  cps  and  At  « 1 , then 

max  3OO 

7500  samples  would  be  required.  A summary  of  these  calculations 
is  shown  in  Figure  3 for  both  the  general  case  and  a specific  example. 


DESCRIPTION  OF  THE  PROGRAM 

The  program  is  written  to  accept  both  gapless  (Datrac)  and  binary 
coded  decimal  (BCD)  tape  input.  Sense  switch  2 governs  the  mode  of 
the  input.  A maximum  sample  of  10,000  words  can  be  handled.  A 
variable  format  for  BCD  data  permits  the  input  to  be  in  any  form 
that  is  specified  on  the  input  format  card.  The  mean  and  standard 
deviation  are  first  calculated  according  to  the  following  standard 
equations: 

mean  = X BAR  1 = - X,  ( 12 ) 


and  printed  out.  It  should  be  noted  that  the  program  does  not 
normalize  the  data  to  have  a zero  mean.  This  cannot  be  done  because 
of  storage  limitations  and  must  be  done  prior  to  input  to  this 
program.  USL  Program  No.  0444  may  be  used  to  normalize  BCD  data. 

Correlation  coefficients  as  a function  of  lag  are  computed 


Standard  deviation  = SIG  X 1 


1 / N 

* [»  (?,  x? 


— 2 
-NX 
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R(p) 


according  to  one  of  the  following  formulas  as  determined  by  sense 
switch  k: 


R(p)  - 


N-p 

(N-P) 

i=l 

N-p  N-p 

% Z** 

1=1  i=l 

K-p 

/N-p  \ 2 1 

N-p 

(N-P)  52  Xi  “ 

Ex 

• 

(»-p)52 

i=l 

\i=l  / J 

i=l 

(14) 


2 

i.+p 


"3  \21|l/2 

Ex 


i +P 


i=l 


N N N 

" Z xi  *hp  -Z  v Z xi 

1=1 1=1 i=l 

r N , _ w_  v2t  r N N p 


(15) 


1/2 


where  p = 0,  1 ...  m and  m is  the  maxinmm  lag.  In  equation  l4  the 

number  of  data  points  used  to  calculate  the  correlation  coefficients 
must  decrease  as  L increases.  In  order  to  use  equation  15,  the  number 
of  points  N,  used  to  calculate  R(p),  remains  fixed  but  the  sample  size 
read  into  the  computer  must  be  greater  than  or  equal  to  N+L.  The 
difference  between  these  equations  should  be  small  because  N will 
be  large  compared  with  m.  Sense  switch  1 will  give  a printed  output 
of  R(p)  and  also  output  on  a plotter  tape. 


Raw  estimates  of  the  power  spectrum  EL  (p)  are  calculated  by  using 
a subroutine  called  Forer.  These  estimates  are  the  result  of  a Fourier 
cosine  transformation  of  the  function  R(p),  and  are  given  by 


EL(p) 


ra-1 

1+2  Z R(q)  008 

q=l 


+ R(m)  cos  np 


p=0,  1, 


m 


(16) 


Finally,  these  raw  estimates  are  smoothed  using  the  Hamming  coefficients 
to  obtain  the  results: 


S(p)  - .23EL(p-l)  + . 54EL(p ) + .23EL(p+l) 

(17) 

where  S(-l)  = S(l)  and  S(M+l)  = S(M-l) 

Sense  switch  6 makes  it  possible  to  obtain  an  average  correlation 
function  for  a set  of  K digitized  sample  functions  {Xj (t  )|  by  finding 
the  arithmetic  mean  of  the  individual  normalized  correlation  functions. 
Power  spectrum  estimates  are  then  found  using  equations  16  and  17  to 
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operate  on  this  average  function.  This  procedui e results  in  a decrease 
in  the  relative  variation  of  the  estimates  of  the  power  spectrum  of  the 
random  process  x(t)  by  a factor  of  approximately  \ '!?,  as  illustrated 
in  Figure  4.  The  solid  curve  represents  the  estimate  c.l  the  power 
spectral  density  of  a band  of  white  Gaussian  noise  as  obtained  from 
.f  single  sample  function,  whereas  the  dotted  curve  wr:;  obtained  by 
averaging  10  samples  in  the  manner  described  above. 

When  using  sense  switch  6,  both  the  individual  c.rreJation  functions 
and  the  average  may  be  printed  out  but  only  the  average  function  arid 
its  transform  are  available  on  a plotter  tape. 

A simplified  flow  chart  given  in  figures  5 and  6 summarizes  the 
order  in  which  all  calculations  are  performed  in  the  program. 


SUMMARY 

An  outline  of  a method  for  estimating  the  power  spectrum  of  a 
stationary  random  px-ocess  from  a digitized  sample  function  as  developed 
by  Blackman  and  Tukey  has  been  given.  A Fortran  program  has  been 
written  utilizing  this  approach.  Parameters  that  must  be  considered 
when  planning  and  before  implementing  such  an  analysis  have  been 
defined  and  emphasis  has  been  placed  upon  the  limitations  imposed  on  the 
resolution  that  can  be  obtained  from  this  type  of  analysis. 
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W(f ) = 2Tra  sin  2nfTm 
2rtfTm 


H(f)  - .5Uw(f)  +.23w/f+_l_\  +.23w/f-  1 ) 

\ 2Tm  / \ 2Tm  / 


Figure  2 
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SAMPLING  RATE 


SAMPLES  / SEC 


RESOLUTION 


DEGREES  OF  FREEDOM 


MAXIMUM  LAG 


NUMBER  OF  LAGS 


APPROXIMATE 
SAMPLE  LENGTH 


APPROXIMATE 
NUMBER  OF  SAMPLES 


lAt=  2fmax 


IDEAL 

LOW 

PASS 


seq 


# = 1/At 


xcps  = I/Tm 


2T ' 
K = Tm 


Tm  = m.At 


M 


T’=l/2K.Tra 


2 ftnax  T' 


EXAMPLE 


W^O  cp^ 


.01  sec 


100/sec 


DESIRE 
.5  cps 


DESIRE 
30  d.f. 


1/.5  = 2sec 


200 


30  sec 


30,000 


PHYSICALLY 
REALIZABLE 
FI LTER 


1 sec 


At-Sfmax 


& m 1/At 


xcps  = I /Tm 


2T* 
K = Tm 


Tm  = M.At  L/.5  = 2sec 


T’=l/2K.Tm 


3 fmax  T* 


EXAMPLE 


fmax=^°  CH 


.0067 


150/sec 


DESIRE 
. 5 cps 


DESIRE 
30  d.f. 


300 


30  sec 


45,000 


Figure  3 
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Figure  4 
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